#include <stdio.h>
#include <math.h>

double habs(double x[],int n)
{
	double max = 0;
	int i;
	for(i = 0;i < n;i++)
		if(fabs(x[i]) > max)
			max = fabs(x[i]);
	return max;
}

void Solve(int n,int disturb)
{
	int i,j,k;
	double h[n][n];
	double b[n];
	double x[n];
	double y[n];
	double t[n];
	//形成 Hilbert矩阵和x向量
	for(i = 0;i < n;i++)
	{
		for(j = 0;j < n;j++)
			h[i][j] = 1.0/(i+j+1.0);
		x[i] = 1.0;
	}
	//形成b向量
	for(i = 0;i < n;i++)
	{
		b[i] = 0;
		for(j = 0;j < n;j++)
			b[i] += h[i][j]*x[j];
		if(disturb == 1)
			b[i] += 0.0000001;
		else if(disturb == 2)
		{
			if(i%2 == 0)
				b[i] -= 0.0000001;
			else
				b[i] -= 0.0000001;
		}
	}
	//形成下三角矩阵
	double l[n][n];
	for(i = 0;i < n;i++)
	{
		for(j = 0;j < n;j++)
			l[i][j] = h[i][j];
	}	
	for(j = 0;j < n;j++)
	{
		for(k = 0;k < j;k++)
			l[j][j] -= l[j][k]*l[j][k];

		l[j][j] = sqrt(l[j][j]);
		for(i = j+1;i < n;i++)
		{
			for(k = 0;k < j;k++)
				l[i][j] -= l[i][k]*l[j][k];
			l[i][j] /= l[j][j];
		}
	}
	//求Lt = b中t值
	for(i = 0;i < n;i++)
	{
		t[i] = b[i];
		for(j = 0;j < i;j++)
			t[i] -= l[i][j]*t[j];
		t[i] /= l[i][i];
	}
	//求L(转置) y = t的y值 
	for(i = n-1;i >= 0;i--)
	{
		y[i] = t[i];
		for(j = n-1;j > i;j--)
			y[i] -= l[j][i]*y[j];
		y[i] /= l[i][i];
	}
	//输出x近似解和x之差的无穷范数
	//printf("x向量和x向量近似解分别为：\n");
	for(i = 0;i < n;i++)
	{
	//	printf("%f %f\n",x[i],y[i]);
		t[i] = x[i] - y[i];
	}
	printf("x之差的无穷范数：\n");
	printf("%f\n",habs(t,n));
	//输出残差向量r和r的无穷范数
	//printf("r残差向量为：\n");
	double r[n];
	for(i = 0;i < n;i++)
	{
		r[i] = 0;
		for(j = 0;j < n;j++)
			r[i] += h[i][j]*y[j];
		r[i] = b[i] - r[i];
		//printf("%f\n",r[i]);
	}
	printf("r的无穷范数：\n");
	printf("%f\n",habs(r,n));		
}

int main()
{
	printf("n = 10，没有扰动：\n");
	Solve(10,0);
	printf("\nn = 10，第一种扰动：\n");
	Solve(10,1);
	printf("\nn = 10，第二种扰动：\n");
	Solve(10,2);
	printf("\nn = 8，没有扰动：\n");
	Solve(8,0);
	printf("\nn = 12，没有扰动：\n");
	Solve(12,0);
	return 0;
}
